#include <string.h>
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "galois.h"
#include "galois_add.h"
#include <map>
#include <sys/time.h>
using namespace std;

#define K 2
#define L 16
#define LK 32
#define Q 65536

int fxn_gcd(uint32_t s, uint32_t j)
{
	uint64_t shift_count, inv, tmpkey, tmpval;
	int i;
	uint64_t Big_Q;
	uint64_t array_size;
	uint64_t high_bit_1, high_bit_2;
	uint32_t high_val_1, high_val_2;
	map <uint64_t, uint32_t> *one;
	map <uint64_t, uint32_t> *two;
	map <uint64_t, uint32_t> *tmpmap;
	map <uint64_t, uint32_t>::iterator it, findit;
	map <uint64_t, uint32_t>::reverse_iterator rit;

	if (LK == 64) Big_Q = 4294967296;
	else Big_Q = Q;

	one = new map <uint64_t, uint32_t> ();
	two = new map <uint64_t, uint32_t> ();

	if (j == 1) array_size = Big_Q+1;
	//else if (j == 2) array_size = Big_Q * Big_Q + 1;

	one->insert(make_pair(1, 1));
	one->insert(make_pair(array_size-1, 1));

	two->insert(make_pair(0, 1));
	two->insert(make_pair(1, s));
	two->insert(make_pair(2, 1));

	rit = one->rbegin();
	high_bit_1 = rit->first;
	high_val_1 = rit->second;
	rit = two->rbegin();
	high_bit_2 = rit->first;
	high_val_2 = rit->second;

	while(1) {
		while(high_bit_1 >= high_bit_2) {
			shift_count = high_bit_1 - high_bit_2;

			inv = galois_composite_divide_k2s(high_val_1, high_val_2, LK, s);

			//for (i = 0; i < array_size; i++) tmp[i] = two[i];
			//for (i = 0; i < array_size; i++) tmp[i] = galois_composite_multiply_k2s(inv, tmp[i], LK, s);
			//for (i = 0; i <= high_bit_2; i++) one[i+shift_count] ^= tmp[i];

			for (it = two->begin(); it != two->end(); it++) {
				tmpkey = it->first;
				tmpval = it->second;

				tmpval = galois_composite_multiply_k2s(inv, tmpval, LK, s);

				findit = one->find(tmpkey+shift_count);
				if (findit == one->end()) {
					one->insert(make_pair(tmpkey+shift_count, tmpval));
				} else {
					findit->second = findit->second ^ tmpval;
					if (findit->second == 0) one->erase (findit);
				}
			}

			rit = one->rbegin();
			if (rit == one->rend()) { break; }
			high_bit_1 = rit->first;
			high_val_1 = rit->second;
		}

		tmpmap = one;
		one = two;
		two = tmpmap;

		if (two->empty()) break;

		rit = one->rbegin();
		high_bit_1 = rit->first;
		high_val_1 = rit->second;
		rit = two->rbegin();
		high_bit_2 = rit->first;
		high_val_2 = rit->second;
	}

	rit = one->rbegin();

	if (rit == one->rend()) return 0;
	high_bit_1 = rit->first;

	if (high_bit_1 == 0) return 1;
	else return 0;
}

main(int argc, char **argv)
{
	uint64_t s, one;
	int i, j, reducible, array_size;
	uint32_t k, z;

	one = 1;

	for (s = 2; s < (one << L); s++) { // f = x^2 + sx + 1
		reducible = 1;
		for (i = 1; i <= K/2; i++) { //only one value for k==2, 2 for k == 4
			/* get (x ^ q ^ i - x) % f */

			k = fxn_gcd(s, i);

			if (k == 0) { reducible = 0; break; }
		}

		if (reducible) {
			printf("%lu\n", s);

			for (j = 0; j < 10000; j++) {
				k = lrand48() % (one << LK);
				z = galois_composite_inverse_k2s(k, LK, s);
				z = galois_composite_inverse_k2s(z, LK, s);

				if (z != k) printf("err with %lu\n", s);
			}
		}
	}
}
